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Abstract 

A Monte Carlo program is presented for the computation of the most gen- 
eral cross section for two-photon production in e + e~ collisions at fixed two- 
photon invariant mass W. Functions implemented for the five 7*7* struc- 
ture functions include three models of the total hadronic cross section and 
the lepton-pair production cross section. Prospects of a structure-function 
determination through a study of the azimuthal dependence between the 
two scattering planes are outlined. All dependences on the electron mass 
and the photon virtualities Q 2 are fully kept. Special emphasis is put on 
a numerically stable evaluation of all variables over the full Q 2 range from 
Qlxam ~ m 2 c (W/y/sY <C m 2 up to Qf max ~ s. A comparison is made 
with an existing Monte Carlo program for lepton-pair production and an 
equivalent-photon approximation for hadronic cross sections. 
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Licensing provisions: 

Computer for which the program is designed and 
others on which it is operable: 
Operating system under which the program has 
been tested: 

Programming language used: 
Number of lines: 
Keywords: 



Monte Carlo, two-photon, e + e , az- 
imuthal dependence 
VEGAS (included, 229 lines) 
RANLUX § (included, 305 lines) 
HBO OK § and DATIME J§ for 
the test program (367 lines) 



none 



all computers 



UNIX 



FORTRAN 77 
1445 



Subprograms used: 



Nature of physical problem: 

Hadronic two-photon reactions in a new energy domain are becoming accessible with 
LEP2. Unlike purely electroweak processes, hadronic processes contain dominant non- 
perturbative components parametrized by suitable structure functions, which are func- 
tions of the two-photon invariant mass W and the photon virtualities Q\ and Q2- It is 
hence advantageous to have a Monte Carlo program that can generate events at fixed, 
user-defined values of W and, optionally, at fixed values of Qi. Moreover, at least one 
program with an exact treatment of both the kinematics and the dynamics over the whole 
range m 2 3> m 2 (W/ \fsY ^ Q\ ^ s (m is the electron mass and y/s the e + e~ cm. energy) 
is needed, (i) to check the various approximations used in other programs, and (ii) to be 
able to explore additional information on the hadronic physics, e.g. coded in azimuthal 
dependences. 
Method of solution: 

The differential cross section for e + e~ — ► e + e~A at fixed two-photon invariant mass W is 
rewritten in terms of four invariants with the photon virtualities Qi as the two outermost 
integration variables in order to simultaneously cope with antitag and tagged electron 
modes. Due care is taken of numerically stable expressions while keeping all electron-mass 
and Qi dependences. Special care is devoted to the azimuthal dependences of the cross 
section. Cuts on the scattered electrons are to a large extent incorporated analytically 
and suitable mappings introduced to deal with the peaking structure of the differential 
cross section. The event generation yields either weighted events or unweighted ones (i.e. 
equally weighted events with weight 1), the latter based on the hit-or-miss technique. 
Optionally, VEGAS can be invoked to (i) obtain an accurate estimate of the integrated 
cross section and (ii) improve the event generation efficiency through additional variable 
mappings provided by the grid information of VEGAS. The program is set up so that 
additional hadronic (or leptonic) reactions can easily be added. 
Typical running time: 

The integration time depends on the required cross-section accuracy and the applied cuts. 
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For instance, 13 seconds on an IBM RS/6000 yields an accuracy of the VEGAS integration 
of about 0.1% for the antitag mode or of about 0.2% for a typical single-tag mode; within 
the same time the error of the simple Monte Carlo integration is about 0.5% for either 
mode. Event generation with or without VEGAS improvement and for either tag mode 
takes about 4 x 10~ 4 (2 x 10~ 3 ) seconds per event for weighted (unweighted) events. 
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1 Introduction 



Two-photon physics is facing a revival with the advent of LEP2. Measurements of two- 
photon processes in a new domain of 77 cm. energies W are ahead of us ||. Any two- 
photon process is, in general, described || by five non-trivial structure functions (two more 
for polarized initial electrons). Purely QED (or electroweak) processes are fully calculable 
within perturbation theory. Several sophisticated Monte Carlo event generators exist 
fTjIIU to simulate 4-fermion production in e + e~ collisions. Indeed, the differential cross 
section is not explicitly decomposed as an expansion in the five 7*7* structure functions. 
Rather, the full matrix element for the reaction e + e~ — > e + e~£ + £~ is calculated as a 
whole, partly even including QED radiative corrections. Such a procedure is, however, 
not possible for hadronic two-photon reactions since the hadronic behaviour of the photon 
is of non-perturbative origin. The decomposition into the above-mentioned five structure 
functions (and their specification, of course) is hence mandatory for a full description of 
hadronic reactions. 

Monte Carlo event generators for hadronic two-photon processes can be divided into 
two classes. Programs of the first kind [|T^,|n|,|T^,|T^jrj|T^| put the emphasis on the QCD 
part but are (so far) restricted to the scattering of two real photons. The two-photon 
sub-processes are then embedded in an approximate way in the overall reaction of e + e~ 
collisions. A recent discussion of the so-called equivalent-photon approximation can be 



found in [16]. 



The other type of programs P7LITBI treat the kinematics of the vertex e + e~ — > e + e~77 
more exactly, but they contain only simple models of the hadronic physics. Moreover, 
the event generation is done in the variables that are tailored for ee — > ee77, namely the 
energies and angles (or virtualities) of the photons and the azimuthal angle between the 
two lepton-scattering planes in the laboratory system. Hence, both the hadronic energy 
W and the azimuthal angle <fi in the photon c.m.s. (which enters the decomposition of 
the e + e~ cross section into the five hadronic structure functions) are highly non-trivial 
functions of these variables^. 

In the study of hadronic physics one prefers to study events at fixed values of W. Not 
only is W the crucial variable that determines the nature of the hadronic physics, but 
through studies of events at fixed W can 77 collisions be compared with 7p and pp ones 
|B| . Next to W, the virtualities Q\ and Q2 of the two photons determine the hadronic 
physics. At fixed values of W and one of the Q's, Qi say, one obtains the cross section 
of deep-inelastic electron-photon scattering. Varying Q 2 one can investigate the so-called 
target-mass effects, i.e. the influence of non-zero values of Q 2 on the extraction of the 
photon structure function F 2 . Hence it is desirable to have an event generator that keeps 
W fixed and in which Qi and Q 2 are the outermost integration variables so that these 
can be held constant. 

The remaining two non-trivial integration variables, which complete the phase space 
of e + e~ — > e + e~X, should be chosen such that three conditions are fulfilled. First, cuts on 
the scattered electrons are usually imposed in experimental analyses. Hence, the efficiency 
and accuracy of the program is improved if these can be treated explicitly rather than 
incorporated by a simple rejection of those events that fall outside the allowed region. 
Second, the peaking structure of the differential cross section should be reproduced as 



1 The only program that contains the ^-dependences is TWOGAM |17 . However, the expressions 
taken from || are numerically very unstable at small Qi, see the discussion following ([lO]). Moreover, 
itself is not calculated. 
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well as possible in order to reduce the estimated Monte Carlo error and to improve the 
efficiency of the event generation. And third, it should be possible to achieve a numerically 
stable evaluation of all variables needed for a complete event description. These three 
conditions are met to a large extent by the choice of subsystem squared invariant masses 
Si and S2 as integration variables besides Q\ and Q\. In the laboratory frame, Sj are 
related to the photon energies Ui by si/ 2 — m 2 = 2u 2 /iy/s, where m denotes the electron 
mass. 

In the interest of those readers not interested in calculational details, the paper starts 
with a presentation of a few results in section |^. The differential cross section for the 
reaction e + e~ — > e + e~X is rewritten in terms of the four invariants Q 2 and Sj (i = 1,2) 
in section ||] where also models for the cross section (7(7*7* — * X) are described. The 
integration boundaries with Q\ and Q\ as the two outermost integration variables are 
specified in section f|. The derivation of the integration limits is standard p(| but tedious. 
Here the emphasis is put on numerically stable expressions^. To our knowledge, numerical 
stable forms of and <fi are presented here for the first time. All dependences on the 
electron mass and the virtualities of the two photons are kept. The formulas are stable 
over the whole range from Qf min ~ m 2 (W / /y / i) 4 <C m 2 up to Q 2 mSuX ~ s, i.e. the program 
covers smoothly the antitag and tag regions. An equivalent-photon approximation is 
also implemented (section ||). The complete representation of the four- momenta of the 
produced particles in terms of the integration variables is given in section Section |7] 
describes the incorporation of cuts on the scattered electrons. Details of the Monte Carlo 
program GALUGA are given in section §. 



2 A few results 

In order to check GALUGA, we include the production of lepton pairs, for which sev- 
eral well-established Monte Carlo generators |7||||J exist. The five structure functions 
for 7*7* — > £ + £~ as quoted in || have been implemented. For the comparison we have 
modified the two-photon part of the four-fermion program DIAG36 || (i.e. DIAG36 re- 
stricted to the multiperipheral diagrams) in such a way that it can produce events at fixed 
values of W. The agreement is excellent. Two examples are shown in Fig. [I], the first 
corresponding to a no-tag setup and the second to a single tagging mode. 

Next we study the (integrated) total hadronic cross section. Figures |2] and |3|-|5] compare 
different ansatze for the Q 2 behaviour of the various cross sections for transverse and lon- 
gitudinal photons. The results of the two models of generalized-vector-meson-dominance 
type (GVMD flU) and VMDc ([15]), dash-dotted and dotted histo grams, respectively) 
are hardly distinguishable in the no-tag case, but may deviate by more than 20% in a 
single-tag case. In the contrast, the different Q 2 behaviour of a simple p-pole (dashed 
histograms) shows up already in the no-tag mode. Note that this model includes scalar 
photon contributions, but does not possess an 1/Q 2 "continuum" term for transverse pho- 
tons. These differences imply that effects of non-zero Q 2 values must not be neglected for 
a precision measurement of a^iW 2 ). 

During the course of the LEP2 workshop, sophisticated programs to generate the 
full (differential) hadronic final state in two-photon collisions have been developed J2I 



The description of hadronic physics with one (or both) photons off-shell by virtualities 

2 A similar phase-space decomposition with s\ replaced by A = [(s — 2m 2 )(W 2 + Q\ + Q\) — {s\ + 
Q\ — m 2 )(s2 + Q\ — m 2 )]/4 is presented in 0. 
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Qj <C W 2 is still premature. Indeed, existing programs are thus far for real photons and 
hence use, in one way or another, the equivalent-photon approximation (EPA) to embed 
the two-photon reactions in the e + e~ environment. It is hence indispensable to check the 
uncertainties associated with the EPA. Hadronic physics is under much better theoretical 
control for deep-inelastic scattering, i.e. the setup of one almost real photon probed by 
the other that is off-shell by an amount Q 2 of the order of W 2 . Corresponding event 
generators exist PI| but also in this case it is desirable to check the equivalent-photon 



treatment of the probed photon. 

An improved EPA has recently been suggested in |16|]. In essence, the prescription 
consists in neglecting Q 2 w.r.t. W 2 in the kinematics but to keep the full Q 2 dependence in 
the 7*7* structure functions. In addition, non-logarithmic terms proportional to m 2 /Q 2 
in the luminosity functions are kept as well. The study [|16[ shows that this improved EPA 



works rather well for the integrated e + e~ cross section. In Fig. ^| we show that this EPA 
(solid compared to dash-dotted histograms) works well also for differential distributions, 
with the exception of the polar-angle distribution of the hadronic system at large angles, 
where it can, in fact, fail by more than an order of magnitude! (There, of course, the 
cross section is down by several orders.) 

The EPA describes also rather well the dynamics of the scattered electrons in the 
single-tag mode except in the tails of the distributions (Fig. |3|). The same holds for the 
distributions in the photon virtualities, see Fig. |. Sizeable differences do, however, show 
up (Fig. H) in the distributions of the subsystem invariant masses ^Jsl. These then lead 
to the wrong shapes for the energy and momentum distributions of the hadronic system 
shown in Fig. |5]. The EPA should, therefore, not be used for single-tag studies. 

Finally we study the prospects of a determination of additional structure functions 
besides F 2 . One such possibility was outlined in 0], namely the study of the azimuthal 
dependence in the 77 c.m.s. between the plane of the scattered (tagged) electron and the 
plane spanned by the beam axis and the outgoing muon or jet. Here we propose to study 
the azimuthal angle between the two electron scattering planes, again in the 77 c.m.s. 
Although such a study requires a double-tag setup, the event rates need not be small, 
since one can fully integrate out the hadronic system but for its invariant mass W . In 
order to demonstrate the sensitivity of such a measurement we show, as a preparatory 
exercise, the distribution for muon-pair production in Fig. []. Fitting to the functional 
form 

— ^ oc 1 + Ai cos + A 2 cos 2 0, (1) 
d0 

we find 

A x = 0.098 , A 2 = -0.028 . (2) 

Let us emphasize that the selected tagging ranges have in no way been optimized for such 
a study. Nonetheless, given the magnitudes of Ai, a measurement appears feasible. 

All but one f?J event generators for two-photon physics use the azimuthal angle be- 
tween the two scattering planes in the laboratory frame as one of the integration variables. 
In fact, appears as a trivial variable in these programs. None of these up to now pro- 
vides the calculation of 0. An expression for in terms of tj, 0, and two other invariants 



is given in || (see ( (49|) below) and, in principle, is available in TWOGAM [fL7| . However, 
the factor \Jt\t2 appears explicitly in the denominator of cos but not in its numerator. 
Hence, at small values of — ij this factor will be the result of the cancellation of several 
much larger terms, rendering this expression for cos numerically very unstable. (Recall 
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that |tj| m i n ~ m 2 (W/y/s) 4 ' <C m 2 , while the numerator contains terms of order s.) In 
contrast, we use the numerically stable expression given in (|50|) P|. 



An approximation for (ft in terms of (ft is proposed in [22 



2 Qi Q2 (2 s - si - s 2 ) 



cos a pp rox = cos (ft + sin (ft . (3) 

( W 2_ tl _ t2 ) / (s _ Sl)(s _ S2 ) 



Indeed, the correlation between (ft and its approximation is very high in the no-tag case, 
where, however, the dependence on (ft is almost trivial (i.e. flat). Figure |6| exhibits that 
there is still a correlation for a double-tag mode, but formula @ fails to reproduce the 
correct (ft dependence: a fit to (P yields A\ = 0.084 and A 2 = 0.017, quite different from 



3 Notation and cross sections 

Consider the reaction 

e + {p a ) + e"(p 6 ) -> e+(pi) + X(p x ) + e~(p 2 ) (4) 
proceeding through the two-photon process 

7(?i)+7(«2) ^X(Px) ■ (5) 

The cross section for (|4]) depends on six invariants, which we choose to be the e + e~ cm. 
energy y/s, the 77 cm. (or hadronic) energy W, the photon virtualities Q i: and the 
subsystem invariant masses J~sl: 



s = (Pa+Pb) 2 , W 2 =p 2 Xl 

81 = (pi + pxf = (pa + qif , -Ql = ti = qf = (po - Pi) 2 , 

s 2 = (P2 + Pxf = (Pb + qif , -Ql = h = ql = (Pb - P2) 2 ■ (6) 

We find it convenient to introduce also the dependent variables: 

2 ( 1 

U2 — si — m — t 2 , v = - [W —h — h 



1 



Ul = s 2 - m 2 - h , K = — ^/X(W 2 , h, t 2 ) = — yV 2 - h t 2 



where X(x, y,z) = (x — y — z) 2 — 4yz and m denotes the electron mass. Note that K is 
the photon three-momentum in the 77 c.m.s. In terms of these variables the e + e~ cross 
section at fixed values of \fs and r = W 2 / s is given by: 

da[e+e- ^ e+ + e" + X] a 2 KW Jn ^ . 2 ~ 2 . 



dr ~ 2^Q 2 Q 2 p 

where R3 is the phase space for (W). 



dR 3 i:(W 2 ,Qi,Qls u s 2}( ft;s } m 2 ) 



3 This form of could, with only minor modifications, be implemented in @]. 
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The hadronic physics is fully encoded in five structure functions. Three of these 
can be expressed through the cross sections o a b for scalar (a, b = S) and transverse 
photons (a, b = T) (o~st = &ts(qi ^ 12))- The other two structure functions ttt and 
tts correspond to transitions with spin-flip for each of the photons with total helicity 
conservation. Introducing 0, the angle between the scattering planes of the colliding e + 
and e~ in the photon c.m.s., these structure functions enter the cross sections as: 

L — 2p 1 Zp 2 a TT + 2pi p 2 CfTS + Pl Z P2 a ST + P\ P2 a ss 

+ 2|p+-p+-|r TT cos2 0- 8 \pt° pt° \ r TS cos0 . (9) 
The density matrices of the virtual photons in the 77-helicity basis are given by 

(«2 — 4 m 2 



00 = (U2 - vf _ , 
Pl K 2 W 2 

\pr\ = pr-i 

I 1 — , /(J30 TT\ l„+-| - "2 



'1 



(pT + 1) \pT\ = i^r Vpi - 1 ( 10 ) 



with analogous formulas for photon 2. 

A few remarks about the numerical stability of the 0-dependent terms are in order. 
Thus far, these terms are implemented solely in the TWOGAM [I7j event generator, 
using the formulas quoted in ||. Given in || and coded in fl7j are the products X 2 = 
2 \p~i~ p 2 ~ I cos 2 and X\ = % \p±° p 2 °\ cos in terms of invariants. Now, the expressions 
for Xi contain explicit factors of tifa (Xz) and \Jt\t2 [X\) in the denominators but not 
in the numerators. Clearly, the evaluation of Xi becomes unstable for small values of 
\ti\. On the other hand, the factors multiplying cos0 and cos 20 in Xi approach perfectly 
stable expressions in the limit m 2 /W 2 — > and ti/W 2 — > 0: 

, +-. 2 „ \ 2m - 2 

\pi + — 



where Xi = W 2 / Si m Sk/s (i 7^ k). Hence a numerically stable evaluation of guarantees 
a correct evaluation of the 0-dependent terms. 

The structure functions o~ a b and r a b for lepton-pair production are well quoted in the 
literature; the formulas of |J are implemented in the program. Much less is known about 
the structure functions for hadronic processes. Since we are not aware of a model for r a b, 
the current version of the program assumes 

T TT = = T TS . (12) 

The cross sections a a b are uncertain at small values of Qi. Three models for a a b are 
provided, all based upon the assumption 

a ab (W 2 , Q 2 ) = h a (Ql) h b (Q 2 2 ) 0-^W 2 ) , (13) 

which is valid for Q 2 <C W 2 , which is justified in most applications. Note the cross section 
for the scattering of two real photons cr 77 (W 2 ) that enters as a multiplicative factor in 
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The three models are defined as follows. The first one is based upon a parametrization 
of the 7*p cross section calculated in a model of generalized vector-meson dominance 
(GVMD): 

h T (Q 2 )=rP 1 -\Q 2 ) + (l-r)P 2 -\Q 2 ) 

m 



hs(Q 2 )=^\r%P l -\Q 2 ) + (l-r 

#s 1 , Q 2 



Q 



2 



2 \nP 2 {Q 2 )-P^{Q 2 ) 
^(<5") = 1 +32> ( 14 ) 



where we take £ = 1/4, r = 3/4, m\ = 0.54 GeV 2 and m 2 = 1.8 GeV 2 . 

The second model |24| adds a continuum contribution to simple (diagonal, three- 
mesons only) vector-meson dominance (VMDc): 

v=£,,p \ m v + Q 2 J m 2 o + Q 2 

hs(Q 2 )= E *^r r v f^TT^V ■ (15) 
v^, P m v Wv + Q 2 ) 

where r p = 0.65, = 0.08, = 0.05, and r c = 1 — J2v r v- 

Since photon-virtuality effects are often estimated by using a simple p-pole only, we 
include also the model defined by (p-pole): 

\m 2 + Q 2 ) m 2 \m 2 + Q 2 J 

Since the program is meant to be used at fixed W, we take 

a^iW) = 1 . (17) 

Finally we give the relation between the cross section at fixed values of r and Q\ and 
the usual form used in deep-inelastic scattering: 

dcr x 2 s da 



drdt 2 Ql dxdQ 2 2 ' v ' 

where x is the Bjorken-x variable defined by 

x = 9l = Ql (ig) 

2 gi .q 2 W 2 + Q 2 + Q 2 - liyj 

4 Phase space 

The phase space can be expressed in terms of four invariants: 

=i6k/ dt2diidSidS2 7^' (20) 
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where A 4 is the 4x4 symmetric Gram determinant of any four independent vectors 
formed out of p a , pb, pi, px, Vi- The physical region in t 2 , t\, si, s 2 for fixed s satisfies 
A4 < 0. Since A4 is a quadratic polynomial in any of its arguments, the boundary of the 
physical region, A 4 = 0, is a quadratic equation and has two solutions. Picking s 2 as the 
innermost integration variable, the explicit evaluation of A4 yields 

16 A 4 = a si + b s 2 + c = a (s 2 - s 2+ ) (s 2 - s 2 _) , (21) 

where 

a = A(si, t 2 , m 2 ) 

b = -2 s mH x - 2 m 2 s\ + 8t 2 m 4 - 2 m 2 t 2 - 2 s s x W^ 2 + 2 m 2 s W 2 + 2 ij s s 1 + 2 s t 2 a x 

t\s — 2 m 2 t 2 s + 2 ti t 2 s\ 



+ 4 m 2 si W 2 + 4 m 4 si + 2 t x t 2 s - 2 1 2 m 2 *i - 2 t 2 s - 2 m 2 t 2 s + 2 i x t 2 si 



- Am 4 W 2 -2ti s 2 + 2s t 2 W^ 2 - 2 m 6 + 2 m 4 ti 

c = -2 s m 4 Vr 2 - 2 t 2 m 2 si - 2 i x t 2 s 2 + 2 s t x t 2 s x - 2 s tjsi + t 2 s 2 + t 2 s 2 + t 2 2 s 2 

+ m 4 s 2 + mH\ - 6 m 6 ti - 2 m 6 si - 4 m 4 si W 2 + 2 m 4 t 2 s + 2 m 4 t 2 s x + 8 m\ s x 

- 2 s 2 t 2 W 2 - 2 ti s 2 - 2 m 2 *! s 2 + m 8 - 2 m 2 s t 2 s 1 +A m 6 W 2 + mH\ 
+ 4 m 2 ti t 2 s - 2 m 2 ti i 2 si - 6 m% + s 2 W 4 + 6 m 2 s t 2 W 2 -As m 2 W 4 
-2s m 2 t\ + 2st 1 m 4 -2s m 2 t\ + 2t 1 t 2 m A + 2s m 2 s 1 W 2 -Ast 1 t 2 W 2 

-2s h m 2 Sl + Qsh m 2 W 2 + 2 s t x si W 2 . (22) 

A numerical stable form for the s 2 limits is 

-b + y/A 

S 2+ = ^ 

2 a 

82- = — , (23) 
a s 2+ 

where A = b 2 — Aac is given below in a numerically stable form, in (fZBj). 

In order to remove the singularity due to (— A4) -1 / 2 (in the limit m 2 <C Sj, W 72 , 
the s 2 integration degenerates to an integration over the <5-function 5(s 2 — sW 2 /s\)), it 
is advisable to change variable from s 2 to x 4 , < x 4 < 1: 



s 2 = — |— b — VK COS (X4 71") j 
ds 2 Air n 



\/-A 4 a/o Jo 



dx 4 . (24) 



For later use we also need a numerically stable form of the Gram determinant, which 
reads 

16A4 = - Asi " 2(X47r) . (25) 
4a 

The si-integration limits follow from the requirement A > 0. They are most easily 
derived when realizing that the discriminant A is given as the product of two 3x3 
symmetric Gram determinants or, equivalently, the product of two kinematic G functions 

^A = AG 3 G 4 = QAD 3 D^ (26) 
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where 

- 4 L> 3 = -4 A 3 (p a , p b , q 2 ) = G(s, t 2 , si, m 2 , m 2 , m 2 ) = G 3 

-4 D 4 = -4 A 3 (p a , gi, g 2 ) = G(t x , s x , t 2 , m 2 , m 2 , PF 2 ) = G 4 . (27) 

Since any 3x3 Gram determinant A3 satisfies A3 > 0, the physical region is that where 
both G3 and G4 are simultaneously negative. Solving Gi for s x 

G 3 = m 2 (s x - s u+ ) (si - su_) 

= m 2 s 2 — 2 m 4 si — s t 2 Si — 3 m 2 t 2 s + m 6 + t 2 s 2 + 
G4 = t x (s x — s X2+ ) (s x — Si 2 _) 

= -2 t x m 2 s x - t 2 m 2 t x + m% - m 2 W 2 ti + m 2 t\ + i 2 VT 2 ti - t x s x W 2 

- 2 m% W 2 + m 2 ^ 4 + t x s\ + i 2 si - t x t 2 s x (28) 



we find 



Sllzt 



t 2 S + 2 m 4 ± V /X(^ m 2 , m 2 ) A(t 2 , m 2 , m 2 ) 



2 m 2 



* 2 2 l^ 2 ti v/A(ti,t 2 ,W 2 )A(ti,m 2 ,m 2 

Si2± = 7T + m + "r; 7T =•= 



2 2 2 2*i 



t 2 ^(-3m 2 + S + t 2 ) 4 
Sii+sii- = 5 hm 



S i 2+ si 2 - = (W^ 2 - m 2 ) {-m 2 + t 2 ) + ^ {W ] ~ hf . (29) 



Note that Si 2+ < si 2 _. Since G 3 is always negative between its two roots, the range of 
integration over s x is Si 2 _ < s x < s xx+ . Numerically it is more advantageous to calculate 
the limits as 



Slmin = 8X2- = TU 2 + - \ W 2 ~ t X + t 2 + Vl yj ' X(W 2 , t X ,t 2 ) 

2 Us + h-Am 2 ) 

Slmax = Sn + = m H — — . (30) 

1 +P V2 

The dominant behaviour of the s x integration is given by the factor \~ 1 / 2 (s x ,t 2 ,m 2 ), 
see (0). (In the limit t 2 , m 2 < sj, this becomes ds x /s x integration.) This factor can be 
transformed away by the variable transformation from si to 23, < £3 < 1, 

si = X x /2 + m 2 + t 2 + 2 m 2 t 2 /X x 
X 1 = (v + KW)(l + y x ) exp{6 x x 3 ) 

5l = ln (v + KW)(l + y x )(l + y 2 ) ' (31) 

such that 

dsi^ = 4tt5i / dx 3 . (32) 



The physical region in the t x -t 2 plane is defined by the requirement Gi < for all s x 
values between the limits (m + W) 2 < s x < {^fs — m) 2 . Since for the reaction considered 
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here the masses of the particles involved are such that the values t\ = (m a — mi) 2 , 
t 2 = {iTLb — 1712) 2 cannot be reached and t 2 is never larger than zero, the boundary curve 
in the t\-t 2 plane is simply given by S\ 2 - = Sn+. Equivalently, the t\ limits can be found 
by solving G4 = with s\ = Sn+ for ti\ 

ilmin — — n \ , ^lmax = 7 , (33) 



where 



ai 

,2 



a 1 = 2(Q + t 2 + 2m 2 + W 2 ) 
b l = Q 2 -W A + 2W 2 t 2 -t 2 2 -8m 2 t 2 -8m 2 W 2 
Cl = 4m 2 (W 2 - t 2 f 
Ai = b\ - AaxCx 

= (Q + t 2 -W 2 + 4mW)(Q + t 2 -W 2 - 4 m W) 

(Q 2 -2Qt 2 + 2QW 2 + t 2 2 + W i -lQm 2 t 2 -2W 2 t 2 ) 

Q = — ^ jt 2 s — m 2 t 2 — m 2 W 2 + \J\{s, m 2 , m 2 )\(t 2 , m 2 , m 2 )| 

- t 2 - W 2 . (34) 



A(s + t 2 -4m 2 ) , , 



1+PV2 

Finally, the ^-integration limits follow from requiring Ai > 0: 

(t 2 — t 2 i + ) (t 2 — t 2 i-) (t 2 — t 22+ ) (t 2 — t 22 J) t 2 (t 2 — t 23 ) 2 . 

Ai ~~ p w p ' ^ > 

p 1 t 2 r 3 



where 



F 2 



As 



t 2 s/3 y 2 + t 2 s - 2 W 2 m 2 + 4 m 3 W 
4s 



t 2 s(3y 2 + t 2 s -2 W 2 m 2 - 4 m 3 W 
F 3 = - 16 m 4 / { -2 1 2 s 2 I3 y 2 + 4t 2 s/3 y 2 m 2 -t 2 s 2 - t 2 s 2 [3 2 + 4 1 2 sm 2 

+4 s 2 f3 2 m 2 - 4 1 2 m 4 + 16 m 6 } 
2 mW -W 2 + s-Am 2 ±P^j(s-W 2 )(s- W 2 ) 



<21± 



^22± 



t 



-2 mW -W 2 + s-Am 2 ±f3 yj(s - W 2 ) (s - W$) 

2 

(s-2m 2 ) 2 



2.3 



m 2 



W± = W ±2m . (36) 
Equivalently, they are arrived at by solving G 3 = with s\ = (m + W) 2 for t 2 : 

t 2 min = t 22+ = ~ [s - W 2 - 2 m W - 4 m 2 + At 2 ) 
m 2 W 2 W% 

^2max — ^22- — T , (37) 

S t-2 m i n 
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where 

At 2 = (3^f(s-W 2 )(s-W 2 ) . (38) 
The phase space finally becomes 

di? 3 = 7-^ / dt 2 / dti 5i(ti,t 2 ) / dxg / dx 4 . (39) 

ifj^ Jtlmin JO JO 

The dominant ^ behaviour is taken into account through a logarithmic mapping, so that 
we end up with a cross section of the form 



dr 



4 r 

n / d^ 

i=l J 

n /dx, ln^ ln^ In- ^ + ^ w r S • (40) 

t 2mi n *1 m i„ (u + K W) f 1 + «1 ) (1 + Wo) 8 IT 2 8 S K ' 



5 Equivalent-photon approximation 

An approximation is arrived at by neglecting as much as possible the electron-mass and 
U dependences in the kinematics, but keeping the full dependence on W and Qi in the 
hadronic cross sections a a b(W 2 , Ql, Q 2 ,) UTS] : 

da _ r s ds 1 r^b dt 2 dti r s ^ ^ ( sW 2 \ a 2 W 2 



dr Jw 2 si Jt 2 a t 2 Jtia t\ Jw 2 \ si J 1Qtt 2 s 

(n „++ o n ++ rr 4- 9 n ++ n 00 rr 
\ r lapprox r2approx U TT T ^ Flapprox r2approx U TS 

+Plapprox 2 P2approx a ST + Pi approx lapprox ^} • (41) 

The integration limits are given by: 



2 w imax 



*ia = ~ . _ - - (1 - a?*) sin o 
^ = -^-(l-x 4 )sin 2 %^, (42) 

\ — X; 2 



where 2i = s 2 /s and a; 2 = S\js. 

The approximate forms of the photon density matrices read: 



2n ++ - 

r lapprox 


2 

x\ 


|l + (l 


„00 _ 
r lapprox 


4 

x\ 





2 2 m 2 x\ 



(43) 



6 Momenta 



Here we present the particle momenta in the laboratory frame. The particle energies 
follow simply from E{ = (p a + p&) • p$/ y^: 

s + m 2 - s 2 

Ei = = 

2^/i 
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E 2 
Ex 



s + m 



Sl 



Si + S2 — 2 m 2 



(44) 



and the moduli of the three-momenta from P 2 = E 2 — mf . The polar angles 9{ with 
respect to the beam axis could be calculated from pt ■ Pi = Ef, Ei — Pb Pi cos 9i 



cos 6*i 
cos 9 2 
cos9x 



s - s 2 + 2t 1 -3m 2 

s - s 1 + 2t 2 - 3 m 2 
2(3^sP 2 
s 2 - si + 2 (t 2 - ti) 
2f3^sP x 



(45) 



Typically, the polar angles are very small and it is better to calculate them in a numerically 
stable form from 



sin 9i 
sin 9 2 
sin 9x 



2^[D[ 

sPPi 

2VD~ 3 

s(5P 2 

2^Dl 

s/3P x 



(46) 



Equations ([4"5D are then only used to resolve the ambiguity 9i «-> n — 6>j. The quantity _D 3 
is defined in (|27| -f29l) ; D\ is obtained from D% by the interchange t\ and s\ «-> S2- 

The same interchange relates -D2, needed below, with .D4, given in (]2"7|--|2"9|) . Furthermore, 
we have: 

D 5 = D 1 + D 3 + 2D 6 



-(s - 4m 2 ) (VT 2 - ^ - t 2 ) + (si - * 2 - ™ 2 ) («a - h - m 2 ) + h t 2 



m 



(si — m 2 ) (s 2 — ^ 2 ) ■ 



(47) 



The polar angles 4>\ (<p 2 ) between the e + (e~) plane and the hadronic plane and the 
polar angle between the two lepton planes in the e + e~ c.m.s. are again best calculated 
using the numerically more stable form for the sinus function 



Da 



COS0 

sin 
sin (f>i 
sin <p 2 



s/3y^Al 

2Vdtd; 



2V=A 4 



s (3 Px sin 9x Pi sin 9\ 

2^Kl 
s (3 Px sin 9x P2 sin 9 2 
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COS (pi 



COS0 2 



£> 3 + As 



r D{ cos i 



-D3 + D\ + 2 y/JjZ^ cos 1 



(4f 



An expression for the azimuthal angle between the lepton planes in the 77 c.m.s. can 
be deduced from the formulas given in B : 



cos 1 



-2 s + Mi + u 2 - u + 4 m 2 + z/(w 2 ~ ") ( Ml ~ ^)/(-^ 2 



2 pi 



(49) 



Numerically more stable is the following form 



sin < 



KW V^AI 



and 



COS 1 



\fD2Dl 



(50) 



where 



16 £> 7 = 2 VT 2 (si s 2 - sW 2 ) - 2 ti (-ti si + sti + si s 2 + W 2 s x - 2 sW 2 

- 2t 2 (-t 2 s 2 + t 2 s + Sl s 2 -2sW 2 + s 2 W 2 ) + 2tit 2 (-si + 2s + 2W 2 
-2 m 2 (m 2 t 2 - t 2 2 - m 2 W 2 + m 2 t x -2W 4 -t 1 2 + W 2 s 1 + 2 t x t 2 
+ 3ti W 2 - t x si + 3t 2 W 2 - t 2 s 2 -t 2 s x + s 2 W 2 - ti s 2 



S-2 



(51) 



A numerically stable relation between and at —ti, m 2 <C W 2 is provided by 

4 s 2 D 2 = 4 si D 3 + 2t 2 s 2 r 2 cos s 2 — 2 ti t 2 ss 2 (s — si) 

— sti t 2 (—2 ii t 2 — st x + t x s x + 3 1 2 s 2 — t 2 s + 2 r 2 cos s) 

— 4 m 2 sr 2 cos si S2 + O (m 4 s 2 s 2 ./ s, m 2 U s s x s 2 ( 

an analogous expression for .D 4 , and 



(52) 



16 s y -D2-D4 cos = 16 y D1-D3 cos — 4 i 2 s 2 r2 cos — 4 t x s 2 r 2 cos 

+ At x t 2 s (— si + 2 s - s 2 + ti + ^2) + 8 m 2 cos r 2 s x s 2 

2 Tfl 2 to / 9 9 9 

+ -* 2 s + 4 st 2 s 2 - 2 i 2 s 2 - 4 ssi s 2 + 2 s x s 2 

s v 

— 8 r 2 cos ss 2 + s 2 s x + s 2 s 2 + 8 r 2 cos s 2 



2 Tfl 2 t-\ 1 > > 

H -ti s + 4 ii ssi - 2 t x si - 4 ssi s 2 + 2 si' 



s 2 



+ r 2 cos </> s ■ 
+ O (m 4 s 2 Sj/s, m 2 t x t 2 



r 2 cos ssi + s 2 Si + s 2 s 2 



where 



2 / S A . . A . *2 si 
m z | — + 1 2 H 

s / V s s 



I — + *i H 

s / V s s 



(53) 



(54) 
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For m — > and ti/W 2 — > 0, (|52| ) and ( [53] ) lead to the approximate relation (|3|). 
The four-momenta are now given by 



Pa — 


' fZ(~\ n n R\ 
- vs [i, u, u, — p) 






Pb = 


~VS(l,0,0,/3) 






Pi = 


(Ei,— Pi sin^i cos 0i. 


, —Pi sin^i sin 0i, - 


-Pi COS^i) 


P2 = 


(#2, --P2 sin ^ 2 cos 2 . 


, P 2 sin0 2 sin 02, Pi 


COS 62) 


Px = 


(E X ,P X sm6 x ,0,P x 


COS^x) • 





7 Experimental cuts 

If cuts on the angle 62 and the energy E 2 of the scattered electron are applied, the (si 
integration region shrinks as follows (see Fig. |7]): 

2 2 1 / 1 2 -^2 max \ 1 



Snow = min |(m + W) , m + s |1 
si U pp = max I ( a/I - m) 2 , m 2 + s (1 



2 E 2 m i n 



and 
where 



max 

) < t 2 < T 2 ( Sl , 2 

min/ 5 



T 2 {si, 6 2 ) — - (3in 2 - s + si + (3 cos #2 sf\(s, si,m 2 )^ 



2m 2 (si -m 2 ) 2 1/2 2 . 2 2 

= uk u72? 2^ 5~^T - A 7 (s, si, m ) sin — 

s [p A 1/Z (s, si, rrr J + s — si — 6 m z \ 2 

_{^ +s(1 _ l2)W |}. 

The approximate form holds for m 2 <C S\ and a small angle 8 2 and is used in (f42|). 
If, as in our case, t 2 is the outer integration, then its lower limit becomes 

^2min = H1U1 {T2(si upp , ^ maK ), Tallow? #2max)} , 

while the upper limit is more complicated 

^2max T2 (-Siupp, 02 min) "§1 -Slupp 

■^llowj 02 min) ^1 ow 

= ^2 Snow < Si < Si U pp , 



where 



2 2 m a/I 
Si = s + m — 



m 2 + 



X 

s(3 2 sin 2 9 2 



X (l + 2m/v / sX 
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2 m 2 — m v^X (0 2 < tt/2) 

2 1 + (3 2 cos 2 fl 2 
2 m —my/s = [U 2 >7r/2) 



I= — + /? 2 sin 2 ^ 2 . (61) 

s 

The Si-integration range is a rather complicated function of £ 2 and may even consist of 
two separated ranges (Fig. |7|). Moreover, the Si-integration range is affected by t 1 and 
cuts on Ei and Q\. Then it is better to use the Monte Carlo method. In any case, since 
the £j integration are the most singular ones, the most important constraints are taken 
into account through ([59]) and (|60| ) and the analogous formulas for t\. 



8 Details of the program 
8.1 Common blocks 

The user can decide whether to keep t 2 at a fixed, user-defined value or to integrate over 
t 2 , i- e - to calculate (|5p or (IB). In the case of integration over all variables, the user can 
choose between the exact or an approximate treatment (|4l| ) of the kinematics. 

Common /ggLapp/ t2user , iapprx , ivegas , iwaght 

t2user Fixed value of t 2 chosen by user for iapprx = 1 . 



iapprx 



ivegas 



ivegas 



t 2 is kept fixed at the user value; 

approximate kinematics is used, t 2 is integrated over; 

all variables are integrated using exact kinematics. 

VEGAS integration; 

Simple integration. 

Unweighted events, i.e. Weight = 1; 

Weighted events 



Cuts on the scattered leptons are set in 

Common /ggLcut/thlmin , thlmax , Elmin , Elmax , th2min , th2max , E2min , E2max 
thlmin,thlmax 



th2min,th2max 



Elmin, Elmax 
E2min,E2max 



Minimum and maximum scattering angles of scattered e^ 
w.r.t. direction of incident e + . 

Minimum and maximum scattering angles of scattered e" 

w.r.t. direction of incident e~. 

Tighter cuts should be applied to the e~. 

Minimum and maximum energies of scattered e + . 

Minimum and maximum energies of scattered e~. 



Models for the 7*7* cross sections and their parameters are chosen in 
Common /ggLmod/ imodel 



imodel 


= 1 


GVMD model (0) 


imodel 


= 2 


VMDc model §L§ 


imodel 


= 20 


p-pole model (|I6D 


imodel 


= 


p-pole model (|B|) with h s {Q 2 ) = 


imodel 


= 3 


Exact cross section for lepton-pair production. 
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Common /ggLhad/ r,xi,mls,m2s,rrho,romeg,rphi,rc,mrhos, 
t momegs ,mphis ,mzeros 



Parameters for (|14}-|T6|): r, £, mi, m|, r p , r^, r^, r c , m pl m 2 , m^,, Wq. 

The integration variables and the particle momenta are stored in 
Common /ggLvar/ 

&yar (4) ,t2,tl , si , s2 ,E1 ,E2,EX,P1 ,P2,PX, thl , th2 , thX , phi 1 , phi2 , phi , pht 

yar(i) Integration variables for VEGAS. 

t2,tl,sl,s2 Invariants t 2 , t±, Si, s 2 . 

E 1 , E2 , EX Energies E u E 2 , E x . 

P 1 , P2 , PX Three-momenta P u P 2 , P x . 

thl,th2,thX Polar angles 6>i, 9 2 , X . 

phil,phi2,phi,pht Azimuthal angles 0i, <p 2 , <P, <P- 

Common /ggLvec/ mntum(7,5) 

Particle four-momenta mntum(i,k): k = 1...5 for p x , p y , p z , E, sign(p 2 ) x J\p 2 \] 
i = 1. . .7 for incident e + , incident e~, photon from e + , photon from e~, scattered e + , 
scattered e~, hadronic system X. 

Parameter for the simple integration and results of the integration and event generation 
are stored in 

Common /ggLuno/ cross, error, Fmax,Fmin, Weight, npts ,nzero ,ntrial 

cross Estimate of luminosity. 

error Estimate of error on luminosity. 

Fmax Maximum function value, calculated in ggLcrs; 

checked in ggLgen. 
Fmin Minimum function value, calculated in ggLuF . 
Weight Weight if weighted events requested, 
npts Number of function evaluations for simple integration, 
nzero Number of cases where function was put to zero in ggLuF 

because it failed the cuts; 

initialized to zero in ggLcrs, ggLgen. 
ntrail Number of trials necessary in ggLgen to generate an event; 

incremented by each call. 

Parameters for the VEGAS integration are set in 

Common /ggLvgl/ xl(10) ,xu(10) ,acc,ndim,nf call,itmx,nprn 

acc VEGAS accuracy (Default (D): 1(T 4 ). 

ndim Number of integration variables (D: 4). 

nf call Maximum number of function calls per iteration for VEGAS (D: 10 5 ). 

itmx Number of iterations for VEGAS (D: 4). 

nprn Print flag for VEGAS (D: 2). 

Additional common blocks 
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Common /ggLprm/ s,roots,Whad,m,Pi,alem 



s Overall cm. energy square s. 

roots Overall cm. energy yfs (twice the beam energy), 

set by user through call to ggLcrs . 
Whad Hadronic mass W, set by user through call to ggLcrs . 
m Electron mass (D: 511 keV). 

Pi 7T 

alem a em (D: 1/137). 

Common/ggLvg2/XI (50 , 10) , SI , SI2 , SWGT , SCHI , NDO , IT 



Common /ggLerr/ 

& itl , iDl , iD3, iD5 , itX, iph, ipl , ip2, ial , ia2 , ia3, ia4 , iel , ie2 , ipt , is 
Block Data ggLblk 



8.2 Subroutines 

ggLcrs (rs,W) Integrates dcx/dr and finds Fmax; rs = y/s, W = W 
Builds up four-momenta. 

Prints four-momenta and checks momentum sum. 



ggLmom 
ggLprt 
ggLgen(Flag) 



Generates one event; 

Flag=F if a new maximum is found; then it is advisable 
to restart event generation with adjusted maximum. 



8.3 Double-precision functions 

ggLint(W2,m2,Qls,Q2s,sl,s2,phi,s) S as defined in 



ggLuF(xar ,wgt) 

ggLhTT(W2,Qls,Q2s) 

ggLhTS(W2,qis,Q2s) 

ggLhSS(W2,qis,Q2s) 

ggLrTS(W2,qis 5 q2s) 

g gLrTT(W2,qis 5 q2s) 

ggLhT(qs) 

ggLhS(qs) 

ggLgg(W2) 

ggLuG(z) 



defined in 



a TT (W 2 ,QlQ 2 2 

°ts(W 2 ,QIQI 

ass(W\QlQl 

r TS (W 2 ,QlQl) 

t tt (W 2 ,QIQI 

h T (Q 2 ) 

h s (Q 2 ) 



{W 2 ) 



Makes the variable transformation from in (EDI) to 
those used by the simple or VEGAS integration. 
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8.4 Excerpt from the demonstration program 

* Initialize the random number generator RanLux 

Call rLuxGo (3,314159265,0,0) 

* 

* Initialize GALUGA; get luminosity within cuts 

Call ggLcrs(rs,W) 

* 

* Initialize plotting 

Call User(0) 

* 

* Timing: 

Call Timex(timel) 

Call rLuxGo (3, 314159265, 0,0) 

* 

* Event loop 

Do 10 i=l,Nev 

Call ggLgen(Flag) 

If ( .not .Flag) Write(6,*) 'Caution: new maximum' 

* 

* Calculate 4-momenta 

Call ggLmom 

* 

* Display first 3 events 

If(i.le.3) call ggLprt 

* 

* Fill histrograms 

Call User(l) 
10 Continue 

* 

Call Timex(Time2) 

Write(6,300) Nev,Time2-Timel , (Time2-Timel)/real(Nev) , 
& iwaght,ntrial,nzero,Fmax 

* 

* Finalize plotting 

Call User(-Nev) 

* 

300 Format (/,3x, 'time to generate ',18,' events is ',E12.5,/, 
&3x, 'resulting in an average time per event of ',E12.5,/, 
&3x, 'unweighted events requested if 1 : ',18,/, 
&3x,'the number of trials was: ',18,/, 
&3x,'the number of zero f was: ',18,/, 
&3x,'the (new) maximum f value was: ',E12.5) 

* 

Stop 
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Figure 1: Comparison of muon-pair production in GALUGA (dashed histograms) and 
DIAG36 (solid histograms) at ^/s = 130 GeV and W = 10 GeV. Top: distribution in the 
logarithm of the polar angle of the fi + [i~ system; no cuts are applied on the scattered 
electrons. Bottom: distribution in t\ under the cuts: 1.55 < Q\ < 3.67° and 30 GeV< E\. 
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Figure 2: Distributions in Ei, ln#i, Ex, and 6x for the integrated total hadronic cross 
section at y/s = 130 GeV and W = 10 GeV. No cuts on the scattered electrons are 
applied. Histogram line-styles correspond to GVMD model in the EPA (solid), p-pole 
model (dashed), GVMD model (dash-dotted), VMDc model (dotted). 
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Figure 3: Distributions in E 1: E 2 , ln^, and ln# 2 for the integrated total hadronic cross 
section at y/s = 130 GeV and W = 10 GeV. The cuts 9 1 < 1.43°, 1.55 < 9 2 < 3.67°, and 
30 GeV< E 2 have been applied. Histogram line-styles correspond to GVMD model in the 
EPA (solid), p-pole model (dashed), GVMD model (dash-dotted), VMDc model (dotted). 
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Figure 4: Same as Fig. [3], but for the distributions in y/sl, ln(— s/ti), and £2 



25 




Figure 6: At the top, the correlation between ap prox ©, proposed in [22H , and 0; at the 
bottom, the distribution in (solid histogram) and its approximation (dashed histogram) 
for the integrated muon-pair cross section at y/s = 130 GeV and W = 10 GeV. The cuts 
1.55° < 6i, 5 GeV< E±, and 30 GeV< Ei have been applied. Also shown is a fit to dcx/d0 
of the form 1 + A\ cos <p + A2 cos 20. 



27 




si 



Figure 7: Phase space in the variables (t 2 ,si) for yfs = 4 and m — 1. The solid lines 
correspond to #2 = 0, 7r/4, 7r/2, 37r/4, and 7r (from t 2 = to t 2 = —12 at Si = 1). The 
dashed lines are si = §1 and t 2 = t 2 at 9 2 = ir/4. 
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